function dphim = conepolyn(wm)
% Calculation of noncommutativity error using polynomial compensation 
% method. Ref. Qin Yongyuan 'Inertial Navigation' P314.
%
% Prototype: dphim = conepolyn(wm)
% Input: wm - gyro angular increments
% Output: dphim - noncommutativity error compensation vector
%
% See also  conetwospeed, conedrift, scullpolyn, cnscl.

% Copyright(c) 2009-2014, by Gongmin Yan, All rights reserved.
% Northwestern Polytechnical University, Xi An, P.R.China
% 09/03/2014, 16/02/2017
    n = size(wm,1);
    if n==1
        dphim = [0,0,0];
    elseif n==2
        dphim = 2/3*cros(wm(1,:),wm(2,:));
    elseif n==3
        dphim = ...
            33/80*cros(wm(1,:),wm(3,:)) + ...
            57/80*cros(wm(2,:),wm(3,:)-wm(1,:));
    elseif n==4
        dphim = ...
            736/945*(cros(wm(1,:),wm(2,:))+cros(wm(3,:),wm(4,:))) + ...
            334/945*(cros(wm(1,:),wm(3,:))+cros(wm(2,:),wm(4,:))) + ...
            526/945*cros(wm(1,:),wm(4,:)) + ...
            654/945*cros(wm(2,:),wm(3,:));
    elseif n==5
        dphim = ...
            123425/145152*(cros(wm(1,:),wm(2,:))+cros(wm(4,:),wm(5,:))) + ...
             34875/145152*(cros(wm(1,:),wm(3,:))+cros(wm(3,:),wm(5,:))) + ...
             90075/145152*(cros(wm(1,:),wm(4,:))+cros(wm(2,:),wm(5,:))) + ...
             66625/145152* cros(wm(1,:),wm(5,:)) + ...
            103950/145152*(cros(wm(2,:),wm(3,:))+cros(wm(3,:),wm(4,:))) + ...
             55400/145152* cros(wm(2,:),wm(4,:));
    elseif n==6
        dphim = ...
            9.225974023727258e-01*(cros(wm(1,:),wm(2,:))+cros(wm(5,:),wm(6,:))) + ...
            8.639610528915165e-02*(cros(wm(1,:),wm(3,:))+cros(wm(4,:),wm(6,:))) + ...
            7.733225109265687e-01*(cros(wm(1,:),wm(4,:))+cros(wm(3,:),wm(6,:))) + ...
            3.930627701652648e-01*(cros(wm(1,:),wm(5,:))+cros(wm(2,:),wm(6,:))) + ...
            5.317640683291427e-01* cros(wm(1,:),wm(6,:)) + ...
            7.627597403941779e-01*(cros(wm(2,:),wm(3,:))+cros(wm(4,:),wm(5,:))) + ...
            3.400757575106209e-01*(cros(wm(2,:),wm(4,:))+cros(wm(3,:),wm(5,:))) + ...
            5.909848488909383e-01* cros(wm(2,:),wm(5,:)) + ...
            7.071861474024891e-01* cros(wm(3,:),wm(4,:));
    else
        dphim = [0,0,0];
%         error('no suitable compensation in conepolyn');
    end